/*
 * BSD Licence:
 * Copyright (c) 2002 Herbert Sauro [ hsauro@cds.caltech.edu ]
 * Systems Biology Group [ www.sys-bio.org ]
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without 
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice, 
 * this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright 
 * notice, this list of conditions and the following disclaimer in the 
 * documentation and/or other materials provided with the distribution.
 * 3. Neither the name of the <ORGANIZATION> nor the names of its contributors
 * may be used to endorse or promote products derived from this software
 * without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
 * DAMAGE.
 */

using System;


namespace Athena.Functions
{

	// ********************************************************************
	/// <summary>
	/// Perform a LU decomposition on the current matrix. The decomposition 
	/// is returned to the caller. The current matrix is not modified.
	/// Code taken from Paul Selormey's article at CodeProject
	/// DotNetMatrix - Simple Matrix library for .NET
	/// http://www.codeproject.com/csharp/psdotnetmatrix.asp
	/// In turn, Selormey's code was taken from
	/// a port of a public domain Java matrix library JAMA. 
	/// 
	/// Code integrated into Matrix class                  
	///
	/// Usage:
	/// <code>
	/// LU = new LUDecomposition (m);
	/// L = LU.getL();
	/// U = LU.getU();
	/// </code>
	/// </summary>
	/// <returns>LU decomposition object is returned. Use getL() and getU() to 
	/// obtain the individual L and U matrices.
	/// </returns>
	// ******************************************************************** 
	public class LUDecomposition {

		/// <summary>
		/// LU matrix
		/// </summary>
		public  Matrix LU;
		/// <summary>
		/// Pivot Sign
		/// </summary>
		public  int pivotSign;
		/// <summary>
		/// Pivot Vector
		/// </summary>
		public int[] pivotVector;

		/// <summary>
		/// Constructor for the LU decomposition object. Use getL() and getU() to obtain the L and U matrices.
		/// </summary>
		/// <param name="m">Matrix to perfom the LU decompoation on.</param>
		public LUDecomposition (Matrix m) {

			this.LU = new Matrix (m);

			// Use a "left-looking", dot-product, Crout/Doolittle algorithm.
			
			int r = LU.r; // m
			int c = LU.c; // n 
			pivotVector = new int[r];

			for (int i = 0; i < r; i++)
				pivotVector[i] = i;
	
			pivotSign = 1;
			double[] LUrowi = new double[c];
			double[] LUcolj = new double[r];
			
			// Outer loop.
			
			for (int j = 0; j < c; j++) {
				
				// Make a copy of the j-th column to localize references.

				for (int i = 0; i < r; i++)
					LUcolj[i] = LU[i,j];
							
				// Apply previous transformations.
				
				// LUrowi is a reference to a row, changes to LUrowi
				// this cause equivalent changes in this[,]
				for (int i=0; i < r; i++) {
					LUrowi = LU.mxx[i];
					
					// Most of the time is spent in the following dot product.
					
					int kmax = System.Math.Min(i, j);
					double s = 0.0;
					for (int k = 0; k < kmax; k++) {
						s += LUrowi[k] * LUcolj[k];
					}
					
					LUrowi[j] = LUcolj[i] -= s;
				}
				
				// Find pivot and exchange if necessary.
				
				int p = j;
				for (int i = j + 1; i < r; i++) {
					if (System.Math.Abs(LUcolj[i]) > System.Math.Abs(LUcolj[p])) {
						p = i;
					}
				}
				if (p != j) {
					for (int k = 0; k < c; k++) {
						double t = LU[p,k]; LU[p,k] = LU[j,k]; LU[j,k] = t;
					}
					int k2 = pivotVector[p]; pivotVector[p] = pivotVector[j]; pivotVector[j] = k2;
					pivotSign = - pivotSign;
				}
				
				// Compute multipliers.
				if (j < r & LU[j,j] != 0.0) {
					for (int i = j + 1; i < r; i++) {
						LU[i,j] /= LU[j,j];
					}
				}
			}
		}


		// ******************************************************************** 
		/// <summary>
		/// Returns lower triangular factor of the LU decomposition. 
		/// 
		/// Usage:
		/// <code>
		/// LU = new LUDecomposition (m);
		/// L = LU.getL()
		/// </code>
		/// </summary>
		/// <returns>Returns the lower triangular matrix</returns>
		// ******************************************************************** 
		public Matrix getL () {
			int m = LU.r;
			int n = LU.c;

			Matrix L = new Matrix(m,n);

			for (int i = 0; i < m; i++) {
				for (int j = 0; j < n; j++) {
					if (i > j) {
						L[i,j] = LU[i,j];
					} else if (i == j) {
						L[i,j] = 1.0;
					} else {
						L[i,j] = 0.0;
					}
				}
			}
			return L;
		}

		// ******************************************************************** 
		/// <summary>
		/// Returns upper triangular factor of the LU decomposition.
		/// 
		/// Usage:
		/// <code>
		/// LU = new LUDecomposition (m);
		/// U = LU.getU()
		/// </code>
		/// </summary>
		/// <returns>Returns the upper triangular matrix</returns>
		// ******************************************************************** 
		public Matrix getU () {
			int m = LU.r;
			int n = LU.c;
			Matrix U = new Matrix(n,n);
			for (int i = 0; i < n; i++) {
				for (int j = 0; j < n; j++) {
					if (i <= j) {
						U[i,j] = LU[i,j];
					} else {
						U[i,j] = 0.0;
					}
				}
			}
			return U;
		}

		// ******************************************************************** 
		/// <summary>
		/// Compute the determinant
		/// </summary>
		/// <returns>Returns a scalar value equal to the determinant</returns>
		// ******************************************************************** 
		public double det () {

			if (LU.r != LU.c)
				throw new EMatrixSizeError ("Matrix must be square to compute the determinant");

			double result = pivotSign;
			for (int i=0; i<LU.r; i++)
				result = result * LU[i,i];
			return result;
		}


	}

}